use PDL; 
use h38b;
use PDL::IO::Dumper;
use PDL::Minuit;
use PDL::Graphics::PGPLOT::Window; 
#my $win = PDL::Graphics::PGPLOT::Window->new(Dev=>'/xs',AspectRatio=>1/1.62,windowwidth=>3);
#use Time::HiRes qw(time);
use strict; use warnings; use feature 'say';
use Storable;

use constant PI => 4* atan2(1,1);
#use constant DEG => PI/180;
use constant DISTANCIA => 2380; # pc
use constant V_LSR => -12; # kms. VLSR default.
my $N_EVALUATIONS = 0;

my $pvdata = rfits("H38b_pv.fits");
my $m1data = rfits("H38b_m1.fits");
$m1data->hdr->{CRVAL1} = $m1data->hdr->{CRVAL2} = 0; # center m1 image to compare with model.
$m1data->hdr->{CTYPE1} = $m1data->hdr->{CTYPE2} = 'OFFSET';
@_=$m1data->dims; ($m1data->hdr->{CRPIX1},$m1data->hdr->{CRPIX2}) = ($_[0]/2,$_[1]/2);
#$m1data->sethdr($m1data->hdr);$m1data->wfits('auxiliar.fits');die;
my $spdata = rfits("H38b_sp.fits");
my $hspd = $spdata->hdr;
my $Velocities = (sequence($$hspd{NAXIS1})+1-$$hspd{CRPIX1})*$$hspd{CDELT1}+$$hspd{CRVAL1};
my $pv_noise = (0.00051505+((yvals($pvdata)>222) | (yvals($pvdata)<146) | 
							(xvals($pvdata)<31) | (xvals($pvdata)>80) )*1e7);
my ($v1mask,$v2mask) = (38e3,56e3);
my $sp_noise = 0.0014+1e7*(($Velocities > $v1mask) & ($Velocities< $v2mask));


# ############## MINIMIZACION ##############
# #_________________________________________
# # Initialization
my @names    =    ('Vlsr','Irot', 'lEM', 'Mass','Incli', # 1 a 5
					'Te3',  'Rd', 'Ri', 'Re', 'Vexp', # 6 a 10
					'offm1x','offm1y','offpvf','offpvv','H'); # 11 a 15
my $pars_ini = pdl(@ARGV);

my ($final_coefs,$final_uncerts,$chi2);
($final_coefs,$final_uncerts,$chi2) = run_minuit({'Run Minimization'=> 0.1,
						  'Minuit Strategy'=>pdl(1),
						  'Fixed Parameters'=>pdl(2,4,6,7,8,11,12,13,14),
						  'Upper Bounds'=>pdl(0,130,0, 0,0,0,0,0,0,0,0,0,0,0,0,0),
						  'Lower Bounds'=>pdl(0, 80,0, 0,0,0,0,0,0,0,0,0,0,0,0,0),
						  'Data'=>{'m1'=>$m1data,'pv'=>$pvdata,'sp'=>$spdata}, 
						  'Noise'=>{'m1'=>3.8,'pv'=>$pv_noise,'sp'=>$sp_noise},
						  #'Cellsizes'=>{'spatial'=>2.5*(DISTANCIA/2380),'velocity'=>1.25} #{au, kms}
						 });
say("Parametros finales ",$final_coefs," -> ",$chi2);
#exit;
plotShimasu($final_coefs,{	'Cellsizes'	=>{'spatial'=>1.25*(DISTANCIA/2380),'velocity'=>1.},
							'Span'	=>{'spatial'=>360*(DISTANCIA/2380),'velocity'=>220},
							'Device'=>'/cps'}); #{au, kms}});


if(0){ # make continuum image
	my ($Vlsr, $ImageRotation, $lEM, $Mass, $Inclination, $Te3, $Rd, $Ri, $Re, $Vexp, 
		    $offm1x, $offm1y, $offpvx, $offpvv, $H) = list($final_coefs);
	my $cube_hash = h38b::cubo({'H'=>$H, 
		    'EM'=>10**$lEM, 'Mass'=>$Mass, 'Incli'=>$Inclination, 'Te'=>$Te3*1e3,  
		    'Rd'=>$Rd, 'Ri'=>$Ri, 'Re'=>$Re, 'Vexp'=>$Vexp,# 
		    'Image rotation'=>$ImageRotation, 'Convolve'=>1, 'Bmin'=>0.020, 'Bmaj'=>0.024, 'Bpa' => 82,
		    'Vlsr'=>$Vlsr,'Distance'=>DISTANCIA,
		    'Cellsizes'=>{'spatial'=>2.5*(DISTANCIA/2380),'velocity'=>1.25} #{au, kms}
		    });
	my $cont = $cube_hash->{'Continuum'}; 
	$cont->sethdr($cube_hash->{'Header continuum'});
	$cont->hdr->{'CTYPE1'} = 'RA---SIN';
	$cont->hdr->{'CTYPE2'} = 'DEC--SIN';
	#print sdump($cont->hdr);
	$cont->wfits("auxiliar.fits");
	system "sethead auxiliar.fits BUNIT=\"Jy\/beam\" ";
}

   
######################################################################
######################################################################
######################################################################

sub run_minuit{
    my ($arg_ref) = @_;
    my $MINIM = $arg_ref->{'Run Minimization'}; defined $MINIM or die;
    my $ubounds = $arg_ref->{'Upper Bounds'} // zeroes($pars_ini);
    my $lbounds = $arg_ref->{'Lower Bounds'} // zeroes($pars_ini);
    my $noise = $arg_ref->{'Noise'}; defined $noise or die;
    my $strategy = $arg_ref->{'Minuit Strategy'} // pdl(1);
    my $verb = $arg_ref->{'Minuit Verbosity'} // pdl(1);
    my $m1data = $arg_ref->{'Data'}->{'m1'};
    my $pvdata = $arg_ref->{'Data'}->{'pv'};
    my $spdata = $arg_ref->{'Data'}->{'sp'};
    my $m1noise = $arg_ref->{'Noise'}->{'m1'};
    my $pvnoise = $arg_ref->{'Noise'}->{'pv'};
    my $spnoise = $arg_ref->{'Noise'}->{'sp'};
    my $fix = $arg_ref->{'Fixed Parameters'};
    my $Cellsizes = $arg_ref->{'Cellsizes'} // {'spatial'=>2.5*(DISTANCIA/2380),'velocity'=>1.25};
    
    my $xi2 = sub{
		my ($npar,$grad,$fval,$xval,$iflag) = @_;
		my ($Vlsr, $ImageRotation, $lEM, $Mass, $Inclination, $Te3, $Rd, $Ri, $Re, $Vexp, 
		    $offm1x, $offm1y, $offpvx, $offpvv,$H) = list($xval);
		my $cube_hash = h38b::cubo_nocont({'H'=>abs($H), 
		    'EM'=>10**$lEM, 'Mass'=>$Mass, 'Incli'=>$Inclination, 'Te'=>$Te3*1e3,  
		    'Rd'=>$Rd, 'Ri'=>$Ri, 'Re'=>$Re, 'Vexp'=>$Vexp,# 
		    'Image rotation'=>$ImageRotation, 'Convolve'=>1, 'Bmin'=>0.020, 'Bmaj'=>0.024, 'Bpa' => 82,
		    'Vlsr'=>$Vlsr,'Distance'=>DISTANCIA,
		    'Cellsizes'=>$Cellsizes}); 
		my $m0model = $cube_hash->{'Moment 0'};$m0model->sethdr($cube_hash->{'Header moment 0'});
		my $m1model = $cube_hash->{'Moment 1'};$m1model->sethdr($cube_hash->{'Header moment 1'});
		my $pvmodel = $cube_hash->{'PV'};$pvmodel->sethdr($cube_hash->{'Header PV'});
		my $spmodel = $cube_hash->{'Flux'};$spmodel->sethdr($cube_hash->{'Header flux'});
		
		my ($m0dr,$m0mr,$m1dr,$m1mr,$pvdr,$pvmr,$arg,$pvnr,$spdr,$spmr,$spnr);
		($m0dr,$m0mr) = h38b::resample($m1data,$m0model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>0});
		return (pdl(1e7),$grad) unless defined $m0mr;
		($m1dr,$m1mr) = h38b::resample($m1data,$m1model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>0});#say "--m1--";
		return (pdl(1e7),$grad) unless defined $m1mr;
		($pvdr,$pvmr,$arg) = h38b::resample($pvdata,$pvmodel,{'Offset'=>pdl($offpvx,$offpvv)});#say "--pv--"; 
		return (pdl(1e7),$grad) unless defined $pvmr;
		$pvnr = $pvnoise->indexND(($arg->{'Pixels first image'})-1);
		#$_ = $pvnr->slice(":,114:125");$_ .= pdl('inf');
		($spdr,$spmr,$arg) = h38b::resample($spdata,$spmodel,{'Offset'=>pdl(0),'Make header'=>1});
		return (pdl(1e7),$grad) unless defined $spmr;
		$spnr = $spnoise->indexND(($arg->{'Pixels first image'})-1);
		my $h = $spdr->hdr; #$v = (sequence($$h{NAXIS1})+1-$$h{CRPIX1})*$$h{CDELT1}+$$h{CRVAL1};#die $v;
		#$m1dr = $m1dr->setbadtoval(0);$m1mr = $m1mr->setbadtoval(0);
		#$m1dr->wfits("auxiliar.fits");die;
		$_=$m0mr->setbadtoval(0);
		my $boolbad = 0;#($m0mr < maximum($_->clump(-1))/10);#die $boolbad->wfits('auxiliar.fits');
		my $chi2s = pdl(sum(($spdr-$spmr)**2/$spnr**2),#->where(($v > 55e3) | ($v < 38e3))), 
				#sum((($m1dr-$m1mr->setbadif($m1dr->isbad))**2/$m1noise**2)->where($m1mr->isgood & $m1dr->isgood)), 
				sum((($m1dr-$m1mr->setbadif($boolbad))**2/$m1noise**2)->where($m1mr->isgood & $m1dr->isgood & !$boolbad)), 
				sum(($pvdr-$pvmr)**2/$pvnr**2));
		my ($NS,$NM1,$NPV) = (sum($spnr<1e7), sum($m1mr->isgood & $m1dr->isgood & !$boolbad), sum($pvnr<1e6));
		#(($pvdr-$pvmr)**2/$pvnr**2)->wfits('auxiliar.fits');die;
		#(($m1dr-$m1mr)**2)->wfits('auxiliar.fits');die;
		#$pvdr->wfits("pvdr.fits");$pvnr->wfits("pvnr.fits");
		#(($pvdr-$pvmr)/$pvnr)->wfits("auxiliar.fits");
		#$chi2s->slice(0) .= $chi2s->slice(0)/nelem($spdr)*4.5;
		#$chi2s->slice(1) .= $chi2s->slice(1)/sum($m1mr->isgood & $m1dr->isgood & $boolbad);
		#$chi2s->slice(2) .= $chi2s->slice(2)/sum($pvmr->isgood)*3;
		$chi2s*=pdl(1/($NS), 1/($NM1+1), 2/($NPV))*($NS/2+$NM1/350+$NPV/2/20)/6;
		$fval = sum($chi2s);
		unless($N_EVALUATIONS%10){
		    #say "$Vlsr, $ImageRotation, $lEM, $Mass, $Inclination, $Te3, $Rd, $Ri, $Re, $Vexp";
		    say "fval: $chi2s -> $fval";
		    #$win->env(-120e3,90e3,-2e-3,30e-3,{xtitle=>'m/s',ytitle=>'Total flux (Jy)'});
		    #$win->points($Velocities,$spdata);$win->line(($cube_hash->{'Velocities'})* 1e3,$spmodel);
		    #$win->line($Velocities->where($spnoise < 100),zeroes($v->where($v < 38e3)),{color=>'black',linewidth=>5});
		    #$win->line($Velocities->where($v > 55e3),zeroes($v->where($v > 55e3)),{color=>'black',linewidth=>5});
		    #$win->line(($cube_hash->{'Velocities'})* 1e3,$spmodel);
		}
		$N_EVALUATIONS++;
		#die;
		return ($fval,$grad);
    };
    
    mn_init(\&$xi2);
    mn_excm('set pri',$verb);# verbosity level
    my $flag = mn_def_pars($pars_ini,.05* ones($pars_ini),{Names=>\@names, Lower_bounds => $lbounds, 
							   Upper_bounds => $ubounds});die "def_pars!\n" if $flag;
    unless($pars_ini->slice(14)){$fix = $fix->append(15)}

    if(defined $fix){$flag=mn_excm('fix',$fix);die "fix!\n" if $flag}
    $flag = mn_excm('set strategy',$strategy); die "set strategy migrad!\n" if $flag;
    $flag = mn_excm('migrad') if $MINIM>=1; die "migrad!\n" if $flag;
    $flag = mn_excm('hesse') if $MINIM; die "hesse!\n" if $flag;
	#$flag = mn_excm('minos') if $MINIM;die "minos!\n" if $flag;
    
    my $chi2final=pdl(mn_stat(1))->slice("(0)");
    @_=mn_pout(1);
    my $final_coefs=pdl($_[0]);
    @_=mn_err(1);
    my $uncerts=pdl($_[2]);
    foreach(1..$#names){
        @_=mn_pout($_+1);$final_coefs=$final_coefs->append($_[0]);
        @_=mn_err($_+1);$uncerts=$uncerts->append($_[2]);
    }
    
    #### Re-normalization  of errors. For each line separatedly line should be ...
    #$uncerts*=sqrt($chi2final/(sum($MODEL_WEIGHTS>0)-nelem($pars_ini)+nelem($fix)));die;
    ####
    unless($MINIM){
        ($chi2final,$_)=\&$xi2(1,zeroes($pars_ini),0,$pars_ini,4);$chi2final=$$chi2final;
    }
    return ($final_coefs,$uncerts,$chi2final);
}

sub plotShimasu{
    my ($final_coefs,$opt_args) = (@_);
    my ($Vlsr, $ImageRotation, $lEM, $Mass, $Inclination, $Te3, $Rd, $Ri, $Re, $Vexp, 
	$offm1x, $offm1y, $offpvx, $offpvv, $H) = list($final_coefs);
    my ($cube,$cube_noconv);
    my $arg_common ={'EM'=>10**$lEM, 'Mass'=>$Mass, 'Incli'=>$Inclination, 'Te'=>$Te3*1e3,  
	      'Rd'=>$Rd, 'Ri'=>$Ri, 'Re'=>$Re, 'Vexp'=>$Vexp,# 
	      'Image rotation'=>$ImageRotation, 
	      'Bmin'=>0.020, 'Bmaj'=>0.024, 'Bpa' => 82, 'Vlsr'=>$Vlsr,
	      'H'=>$H,
	      'Distance'=>DISTANCIA,
    };
    $cube = h38b::cubo({ %$arg_common , %$opt_args, 'Convolve'=>1});
    
    my $m0model = $cube->{'Moment 0'};$m0model->sethdr($cube->{'Header moment 0'});
    my $m1model = $cube->{'Moment 1'};$m1model->sethdr($cube->{'Header moment 1'});#die sdump($cube->{'Header moment 1'});
    my $pvmodel = $cube->{'PV'};$pvmodel->sethdr($cube->{'Header PV'});
    my $spmodel = $cube->{'Flux'};$spmodel->sethdr($cube->{'Header flux'});
    
    #die $spmodel;
    #$cube->{'Header PV'}->{VRESOL} = 1.35; # for the error bar resolution in the pv diagram
    #$pvdata->hdr->{VRESOL} = 1.35;
    
    my ($m0dr,$m0mr,$m1dr,$m1mr,$pvdr,$pvmr,$arg,$pvnr,$spdr,$spmr,$spnr);
    ($spdr,$spmr,$arg) = h38b::resample($spdata,$spmodel,{'Offset'=>pdl(0),'Make header'=>1});
    $spnr = $sp_noise->indexND(($arg->{'Pixels first image'})-1);
    #die $spnr-$spdr;
    my $h = $spdr->hdr; #say sdump($h);
    ($m1dr,$m1mr) = h38b::resample($m1data,$m1model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>1});
    ($m0dr,$m0mr) = h38b::resample($m1data,$m0model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>1});
    #say sdump($m1dr->hdr);$m1dr->wfits("auxiliar.fits");
    ($pvdr,$pvmr,$arg) = h38b::resample($pvdata,$pvmodel,{'Offset'=>pdl($offpvx,$offpvv),'Make header'=>1});
    $pvnr = $pv_noise->indexND(($arg->{'Pixels first image'})-1);

    my $vels = (sequence($$h{NAXIS1})+1-$$h{CRPIX1})*$$h{CDELT1} + $$h{CRVAL1};#die (nelem($vels),nelem($spdr));
    #say sdump($m1dr->hdr);die;
    $pvdr->hdr->{VRESOL} = 1.35; # for the error bar resolution in pv
    $pvdr->hdr->{B} = 5.56e-6* 3600; # for the error bar resolution in pv
    
    my $max_pvdr = maximum($pvdr->clump(-1));#die $spmr;
    use PDL::ImageND;

    my $device = $opt_args->{'Device'} // "/xs";
    if($device =~ /ps/){$device = "plota.ps$device"}
    my ($vl1,$vl2)=(-37,13);
    h38b::triplePlotData({'Velocities'=>$vels/1e3,'Flux'=>[$spmr,$spdr],
			  'Moment 1'=>$m1dr,'Header moment 1'=>($m1dr->hdr),
			  'Moment 0'=>$m0dr,'Header moment 0'=>($m0dr->hdr),
			  'PV'=>convolve($pvdr,ones(3,3)/9),'Header PV'=>($pvdr->hdr),
			  'Image rotation'=>$ImageRotation},
			 {'Draw beams'=>1,'Layout'=>'horizontal',
			  'PV contours'=>(sequence(4)+1)*.2*$max_pvdr,
			  'M1 limits'=>pdl($vl1,$vl2),'M1 contours'=>$vl1+($vl2-$vl1)*sequence(6)/5,
			  'M0 no contours'=>1, 'Device'=>$device
			  ,'PV Label'=>["Data",.13,.80]
			 });
    print "PV contours max: $max_pvdr\n";
    $h = $m1mr->hdr;
    #$m1mr = $m1mr->setbadif($m1dr->isbad);$m1mr->sethdr($h);
    $_ = $m0mr->setbadtoval(0); my $boolbad = ($m0mr < maximum($_->clump(-1))/10);
    $m1mr = $m1mr->setbadif($boolbad);$m1mr->sethdr($h);
    if($device =~ /ps/){$_ = $device;s/plota/plotb/;$device = $_}
    h38b::triplePlotData({'Velocities'=>$vels/1e3,'Flux'=>$spmr,
			  'Moment 1'=>$m1mr,'Header moment 1'=>($m1mr->hdr),
			  'Moment 0'=>$m0mr,'Header moment 0'=>($m0mr->hdr),
			  'PV'=>$pvmr,'Header PV'=>($pvmr->hdr),
			  'Image rotation'=>$ImageRotation},
			 {'Draw beams'=>1,'Layout'=>'horizontal',
			  'PV contours'=>(sequence(4)+1)*.2*$max_pvdr,
			  'M1 limits'=>pdl($vl1,$vl2),'M1 contours'=>$vl1+($vl2-$vl1)*sequence(6)/5,
			  'M0 no contours'=>1, 'Device'=>$device,'PV Label'=>["Model",.13,.80]
			 #'M1 cutoff'=>4.5*0.000741922692 # in Jy/bm
			 });
	if($device =~ /ps/){$_ = $device;s/plotb/plotd/;$device = $_}
	#$m1dr->wfits('auxiliar.fits');die; 
    h38b::triplePlotData({'Velocities'=>$vels/1e3,'Flux'=>$spdr-$spmr,#'Flux weight'=> $spnr<1e6,
			'Moment 1'=> ($m1dr-$m1mr)->setbadif($m1mr->isbad | $m1dr->isbad),'Header moment 1'=>($m1mr->hdr),
			'Moment 0'=>$m0dr,'Header moment 0'=>($m0mr->hdr),
			'PV'=>(($pvdr-$pvmr)->setbadif($pvnr>=1e6)),'Header PV'=>($pvmr->hdr),
			  'Image rotation'=>$ImageRotation},
			 {	'Flux limits' =>[-7e-3,7e-3],
			 'Flux mask'=>{'Limites'=> pdl($v1mask/1e3,$v2mask/1e3,-7e-3,7e-3),'Filling'=>'hatched','Color'=>'black'},
			 	'Draw beams'=>1,'Layout'=>'horizontal',
			  'PV contours'=>"colormap",
			  'M1 limits'=>pdl(-15,15),
			  'M0 no contours'=>1, 'M1 no contours'=>1, 
			  'Device'=>$device, 'M1 colormap'=>'aips0',
			  'PV Label'=>["Data - Model",.1,.85,{charsize=>2,colour=>'white'}]
			 #'M1 cutoff'=>4.5*0.000741922692 # in Jy/bm
			 });

    $$arg_common{'Re'} *=1.0; 
    
    #$cube_noconv = h38b::cubo({ %$arg_common , %$opt_args, 'Convolve'=>0, Cellsizes=>{spatial=>1.5,velocity=>3}});
    $cube_noconv = h38b::cubo({ %$arg_common , %$opt_args, 'Convolve'=>0});
    $m0model = $cube_noconv->{'Moment 0'};$m0model->sethdr($cube_noconv->{'Header moment 0'});
    $m1model = $cube_noconv->{'Moment 1'};$m1model->sethdr($cube_noconv->{'Header moment 1'});$m1model->wfits('m1model.fits');
    #$m1model->wfits('auxiliar.fits');die;
    $pvmodel = $cube_noconv->{'PV'};$pvmodel->sethdr($cube_noconv->{'Header PV'});
    $spmodel = $cube_noconv->{'Flux'};$spmodel->sethdr($cube_noconv->{'Header flux'});
   
    ($m0dr,$m0mr) = h38b::resample($m1data,$m0model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>1});
    ($m1dr,$m1mr) = h38b::resample($m1data,$m1model,{'Offset'=>pdl($offm1x,$offm1y)/3600,'Make header'=>1});
    ($pvdr,$pvmr,$arg) = h38b::resample($pvdata,$pvmodel,{'Offset'=>pdl($offpvx,$offpvv),'Make header'=>1});
    ($spdr,$spmr,$arg) = h38b::resample($spdata,$spmodel,{'Offset'=>pdl(0),'Make header'=>1});
    if($device =~ /ps/){$_ = $device;s/plotd/plotc/;$device = $_}
    my $max_pvmodel = $pvmodel->clump(-1)->maximum;
    h38b::triplePlotData({'Velocities'=>$vels/1e3,'Flux'=>$spmr,
			  'Moment 1'=>$m1model->setbadtoval(-pdl('inf')),'Header moment 1'=>($m1model->hdr),
			  #'Moment 0'=>$m0mr,'Header moment 0'=>($m0mr->hdr),
			  'PV'=>convolveND($pvmodel,rvals(3,3)/sum(rvals(3,3))),'Header PV'=>($pvmodel->hdr),
			  'Image rotation'=>$ImageRotation},
			 {'Draw beams'=>0,'Layout'=>'horizontal',
			 'M1 no contours'=>1,
			 #'M1 contours'=>pdl($Vlsr-sqrt(2)*3.59*$Te3**.5,$Vlsr+sqrt(2)*3.59*$Te3**.5), # \mu_m = 0.64 m_H
			'M0 no contours'=>1,
			 'PV contours'=>(sequence(4)+1)*0.2*$max_pvmodel,
			 'M1 limits'=>pdl(-35,15), 'Device'=>$device,'PV Label'=>["Model",.13,.80]
			 });
    if($device =~ m/ps$/){
    	system "ps2eps -B -f plota.ps; epstopdf plota.eps; pdf270 plota.pdf -o plota.pdf";
    	system "ps2eps -B -f plotb.ps; epstopdf plotb.eps; pdf270 plotb.pdf -o plotb.pdf";
    	system "ps2eps -B -f plotc.ps; epstopdf plotc.eps; pdf270 plotc.pdf -o plotc.pdf";
		system "ps2eps -B -f plotd.ps; epstopdf plotd.eps; pdf270 plotd.pdf -o plotd.pdf";
    	system "rm -f plot*ps";
    }
}

